rm(list = ls())

# Functions
R <- function(x){
  out <- pbeta(q = x, shape1 = 2, shape2 = 2) - 
    dbeta(x = x, shape1 = 2, shape2 = 2)*x
  return(out)
}

x.0 <- function(){
  out <- uniroot(R, lower = 0.5, upper = 1)$root
  return(out)
}

# Figure
labels.x1 <- c('0', '0.25', expression(hat(x)), 0.75, '1')
labels.x2 <- c('0', '0.25', expression(hat(x)), expression(x[0]), '1')
at.x1 <- c(0, 0.25, 0.5, 0.75, 1)
at.x2 <- c(0, 0.25, 0.5, x.0(), 1)

X <- seq(0, 1, length.out = 1000)

pdf('./figures/fig_c1.pdf', height = 7, width = 14)
par(mfrow = c(1, 2))

plot(X, dbeta(X, shape1 = 2, shape2 = 2), 
     type = 'l', lwd = 5, xaxt = 'n', 
     ylab = 'Probability / Density', 
     xlab = 'x', 
     cex.lab = 1.35, cex.axis = 1.35)
axis(side = 1, labels = labels.x1, at = at.x1, 
     cex.axis = 1.35)
lines(X, pbeta(q = X, shape1 = 2, shape2 = 2), 
      col = 'gray', lwd = 5)
segments(x0 = 0.5, y0 = 0, 
         x1 = 0.5, 
         y1 = dbeta(0.5, shape1 = 2, shape2 = 2),
         lwd = 3, lty = 2)

plot(X, R(x = X), 
     type = 'l', lwd = 5, xaxt = 'n', 
     ylab = 'R(x)', xlab = 'x', 
     cex.lab = 1.35, cex.axis = 1.35)
axis(side = 1, labels = labels.x2, at = at.x2, 
     cex.axis = 1.35)
abline(h = 0, lwd = 3, lty = 2)
segments(x0 = 0.5, y0 = 0, 
         x1 = 0.5, 
         y1 = R(0.5),
         lwd = 3, lty = 2)
segments(x0 = x.0(), 
         y0 = -0.25, 
         x1 = x.0(), 
         y1 = R(x.0()),
         lwd = 3, lty = 2)

dev.off()